#include<iostream>
using namespace std;
#include<iomanip>
#include<ctime>
#include<cmath>
#include<omp.h>

#define OMP_NUM_THREADS 10

/*
 * integ_midpoint: 中值法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_midpoint(double a, double b, int n, double (*func)(double)){

	double y=0;
	double h=(b-a)/n;
	#pragma omp parallel for reduction(+:y)
	for(int i=0;i<n;i++){
		y += (a+ (i+0.5)*h )*(a+(i+0.5)*h);
	}
	return y*h;

	/*
	omp_set_num_threads(OMP_NUM_THREADS);

	#pragma omp parallel
	{
	       double x,sum;
	       int i, id, num_threads;
	       id = omp_get_thread_num();
	       num_threads = omp_get_num_threads();
	       
	       for(i=id,sum=0.0; i<n;i= i+ num_threads){
		       x = (a + (i+0.5)*h );
		       sum += func( a + (i+0.5)*h );
	       }
		#pragma omp atomic
	       y += sum;
	}

	return y*h;
	*/
}

double f(double x){
	return x*x;
}

int main(){

	double a=0,b=1, y;
	int n=1E7;

	clock_t t_start, t_end;
	t_start = clock();
	y=integ_midpoint(a,b,n,f);
	t_end = clock();
	cout<<"Elapsed time: "<<setprecision(15)<<(double)(t_end-t_start)/CLOCKS_PER_SEC<<" s.\n";
	cout<<"integ = "<<y<<endl;

	return 0;
}
